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Abstract 

Using Langevin dynamics simulations, we study elastic response of single semiflexible polyelec- 
trolytes to an external force pulling on the chain ends, to mimic the stretching of DNA molecules 
by optical tweezers. The linear chains are condensed by multivalent counterions into toroids. The 
force-extension curve shows a series of sawtooth-like structure, known as the stick-release patterns 
in experiments. We demonstrate that these patterns are a consequence of the loop-by-loop unfold- 
ing of the toroidal structure. Moreover, the dynamics, how the internal structure of chain varies 
under tension, is examined. At the first stage of the stretching, the toroidal condensate decreases 
its size until the loss of the first loop in the toroid and then, oscillates around this size for the 
rest of the unfolding process. The normal vector of the toroid is pulled toward the pulling-force 
direction and swings back to its early direction repeatedly when the toroidal chain looses a loop. 
The results provide new and valuable information concerning the elasticity and the microscopic 
structure and dynamic pathway of salt-condensed DNA molecules being stretched. 

PACS numbers: 82.35.Rs, 87.15.La, 87.15.hp, 87.15.ap 



1 



There are continually strong demands in understanding the properties of DNA molecules 
because of their broad implications in biology and the benefits in gene therapy l|. Owing 
to the progress of nano-technology, researchers are now allowed to investigate the elastic 
response of DNA under the action of an external force, at the level of single molecules . 



When a DNA mo 
been observed 



ecule was stretched, different kinds of force-extension curve (FEC) have 
% 15J, depending on the multivalent counterion concentration, C, in the 
solution. For C staying outside the region bounded by the condensation threshold C c and 
the decondensation threshold Cd (C c < Cd), DNA molecules are in a coil state and FEC is 
basicallv described by the Marko-Siggia formula derived from the wormlike chain (WLC) 
model 6|. For C inside the region, the DNAs are collapsed into globule of ordered structure 
and the preferable structure is toroid [7j. At this moment, FEC shows plateau force while 
C ~ C c , or stick-release pattern while C is close to Co, where Co is the salt concentration 
at which the effective chain charge is neutralized by the multivalent counterions [8|. For the 
plateau force, the elastic response is constant over a wide range of chain extension and its 
value is larger than the prediction of the WLC model, whereas for the stick-release pattern, 
the force is periodic and piecewisely described by the Marko-Siggia formula j^l, 5]. 

It was suggested that the plateau force is derived from the fact that the ratio of the 
extension to the effective chain contour length is a constant and thus produces a constant 
force according to the Marko-Siggia formula M. To explain the stick-release pattern, two 
models have been proposed. The first model [9|, llO] is based upon the observation of the 
rings- on- a-string structure of a long DNA chain. The formation of this gripping pattern 
could be a consequence of the abrupt breaking of one ring into several under the stretching 
of an external force. Elastic response between any two ring breakings is determined by the 
coil part of the chain and follows the WLC model. The second model , on the other 

hand, attributes this phenomenon to a result of the loop-by-loop unfolding of a toroidal 
DNA condensate under stretching, because a pronounced quantization in the DNA release 
length of about 300nm has been demonstrated, which correlates exactly with the periphery 
length of a typical DNA toroid of size R ~ 50nm. Nevertheless, a concrete evidence is still 
missing and researchers cannot make a conclusive connection of the stick-release pattern 
with the second model yet. 

Recently, the pathway of unwrapping a spool of DNA chain helically coiled onto a histone 
protein has been studied by theorists to investigate the stability and dynamics of nucleo- 
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somes under tension [12j. They pointed out the similarity between two seemingly unrelated 
problems: unwrapping of nucleosomes and unfolding of DNA toroidal condensates, and pre- 
dicted a catastrophic event for the two systems under tension: a sudden and quantized 
unraveling happened once a time for a DNA turn. This unwrapping pathway is difficult 
to be investigated by experiments and the information concerning the internal structure 
of toroidal DNA condensates under tension is still not complete. Therefore, in this study, 
we conduct computer simulations to study the elastic response of a toroidal polyelectrolyte 
(PE) condensate and investigate structural variation and unfolding dynamics of chain un- 
der tension. We focus on the case that the amount of multivalent counterions is in charge 
equivalence with the PE chain. This choice gives us the most chance to observe the stick- 
release pattern. There have been simulation works devoting to the study of PE chains under 



tension 
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141 ] . However, in most of the works, the chains were flexible and condensed 



into disordered structures either by compacting potentials or by monovalent counterions 
with a ultra-strong Coulomb coupling. Our work here, on the contrary, deals with a more 
realistic situation in an ambient condition where the chain is semiflexible, collapsed by mul- 
tivalent counterions into a toroid. This kind of study is still few 15] , and involves nonlinear 
and nonequilibrium effects coming from the arrangement of ions, the long-range Coulomb 
interaction, and the pulling speed. 

Our system contains a single chain and many counterions. The chain is constituted of 
N m = 512 monomers, each of which carries a negative unit charge — e. The counterions are 
tetravalent and the number of the counterions is 128. The excluded volume of the monomers 
and the counterions is modeled by the Lennard- Jones potential £lj [2(c/^) 6 — l] 2 truncated 
at the minimum. We assumed an identical e^j for the monomers and counterions but 
different diameter a. The diameter of the counterions a c is half of that of the monomers o m . 
In the following text, we use £lj and a c as the units of energy and length, respectively. Two 
consecutive monomers on the chain are connected by a bond of length b via the harmonic 
potential kf,(b — bo) 2 with kf, = 100 and b$ = 1.1. The chain stiffness is described by an 
angle potential /^(^ — #o) 2 + k^(9 — 9q) a with kz = 5, = 200 and 9q — where 9 is 
the angle between two adjacent bonds on the chain. Coulomb interaction is expressed by 
X^ksTziZj/r where r is the separation distance of two particles of valences Zi and Zj, and 
Ab is the Bjerrum length. We set Ab = 4.68 and the temperature k B T = 1.2. We assumed 
further that the bond potential and the angle potential have already incorporated the effect 
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of non-bonded interactions along the chain, including the excluded volume interaction and 
Coulomb interaction. Therefore, the non-bonded interactions are deactivated for pairs of 
monomers separated by less than three bonds on the chain. This kind of setup has been 
used in simulations. The advantage of this setup is that the bond and the angle potentials 
given here have been the entire potentials without need to take into account further the non- 
bonded interactions between the monomers. Please notice that the deactivation of Coulomb 
interaction locally on the chain will not violate the requirement of charge neutrality because 
it is a global property of a system and the total charge is always zero. The system is 
placed in a rectangular box of size 620 x 91 x 91 and periodic boundary condition is applied. 
The technique of PPPM Ewald sum is employed to calculate the Coulomb interaction. We 
performed Langevin dynamics simulations [16] with the friction coefficient ( setting to 2t~ 1 



and It -1 for the monomers and the counterions, respectively, where r = o c \jmje^ is the 
time unit and m is the particle mass. Stevens has used a similar model [lTj and shown that 
a PE chain can collapse into a toroid due to the condensation of tetravalent counterions. 
Since the line charge density matches that of a double-stranded DNA (dsDNA) and the size 
of monomer matches the size of a phosphate group, our model can be used to understand 
the system of DNA condensed into a toroid. We remark that the persistence length in our 
model is one order of magnitude shorter than dsDNA. The reason for this setup is to allow 
the formation of a toroid in a short chain length; thus the simulations become feasible under 
limited computing resources. 

The initial configuration of chain is a circular helix with its central line lying on the x axis. 
This configuration is used to favor the formation of a toroid. Other structures, such as cigar 
or racket shapes, sometimes appear in the simulations 17lll8l|. All these structures have been 
observed in experiments of DNA condensation and people believe that toroid is the most 

n 

stable structure [7]. Since different condensed structure could produce different FEC under 
tension, in this study we focus on the case of chains collapsing into toroids. The two chain 
ends are guided by a spring force towards two points on the x-axis with separation distance 
equal to 60. After reaching a stable state, most of the chains form rod-toroid-rod structure 
as shown in the inset of Fig. [IJ Since chain ends are constrained, knots will not appear on 
the chain. We have chosen different radii and pitches of helix as our initial configurations 
and verified that the size of the generated toroid is independent of these choices. The toroid 
size is defined by two radii, the minor radius ro and the major radius Rq. The previous one 
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is the radius of the annular tube of the toroid and the latter is the radius measured from 
the toroid center to the center of the annular tube. These two radii determine the gyration 
tensor of a toroid with the three eigenvalues equal to rjj/4, (4i?Q + 3rQ)/8, and (4Rl + 3rl)/8. 
The gyration tensor of a set of N particles can be calculated by simulations according to the 
equation T af3 = X)i=i(^ — r cm ) a {fi — r cm )p/N where r cm is the center of mass of the set of 
particles and the subscripts a and j3 denote the three Cartesian components. We calculated 
the gyration tensor of the toroid on the chain and estimated vq and Ro from the three 
eigenvalues A, (i = 1,2,3). The results are r = 3.60(1) and i? — 12.3(1). It is known that 
a ring structure (r <C Ro) can be characterized by a quantity, called asphericity, defined 
as A = Yli=i Sj=i(Ai — Aj) 2 / (2 Yli=i Aj) 2 , with its value equal to 0.25. Since our toroid has 
a finite ro, A is smaller than 0.25 and takes a value of 0.235(1). The winding number W 
is another important quantity to describe the status of a toroid, which counts the number 
of turns winding around the toroid central axis. Our toroidal condensates have W = 6.5 
before stretching. The decimal 0.5 in W reflects the fact that the two chain ends stay on 
the opposite sides of the toroid. 

We fixed one end of the condensed PE chains and pulled the other end outwards along 
the x-axis using a spring force with constant pulling speed v = 0.0005. Since unfolding of a 



condensed PE is a nonequilibrium process, it is delicate to choose the pulling speed [U|, 
We have verified that significant nonequilibrium effect is detected only when the pulling 
speed excesses the characteristic speed, v s = Rq/tr ~ 0.002, estimated by the Rouse re- 
laxation time tr. Our choice of v is slow enough to minimize the dependence of FEC on 
the pulling speed, which essentially probes the elastic response of PEs in the limit of zero 
pulling speed [l4]. Nevertheless, it does not mean that our chains reach equilibrium for ev- 
ery stretched distance. To study it, one needs to perform equilibrium simulations with the 
chain ends fixed on every stretched distance. However, this sort of approach demands much 
more computing resources. For practical purpose to understand the stick-release pattern 
of a toroidal DNA condensate under tension, we decided to follow the arguments given by 
Lee and Thirumalai [jjj] that the nonequilibrium method with a small pulling speed can 
be a good approach, sufficient enough to study the problems of the elasticity, the internal 
structure, and the dynamics. 

We prepared 20 independent toroids and performed stretching process. A typical FEC 
obtained in our simulations is shown in Fig. [TJ The extension x is the distance between 
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FIG. 1: Force / vs. extension x (thick solid curve). The dashed line is the fitting curve for x > 526 
by the EWLC model. W is plotted in solid curve and the value is read from the right axis of 
ordinate. The inset is a snapshot of the chain before stretching. 



two chain ends. The data have been passed through a low-pass filter to attenuate the high- 
frequency noise. We have verified that the filter does not modify the content of the force 
curve. The FEC shows approximately a plateau force up to x ~ 270. After that, the force 
increases and then decreases abruptly with extension. This tooth-like structure repeatedly 
appears in the FEC and becomes shaper and sharper as the extension increases. This is 
so-called stick-release pattern, which has been experimentally reported 3|, |5[]. We plot in the 
same figure the winding number W. We see that W is a downstairs function and the height 
of each stairstep is equal to one. It shows that the toroidal PE looses one loop, followed 
by one loop, under stretching. The results support the picture of Murayama et al. 5]. The 
dynamics that one loop breaks into several was not observed in this study. Moreover, we 
observed that the sawteeth in FEC appear coherently with the steps of W . Each sawtooth 
corresponds exactly to one step in W. It demonstrates that the stick-release pattern is a 
consequence of the loop-by-loop unfolding of the condensed PE toroid, and confirms the 
theoretical picture of stepwise unwinding of PE under stretching 191 ]. 

In the final stage of the stretching (x > 525), / increases largely and the elastic response 
is similar to the WLC model. Since our chain is extensible, we fitted these data by the 
modified Marko-Siggia formula 20j derived from the extensible WLC (EWLC) model, 



fP 
k B T 



x 1 1 

L + 4(1- x/L + f/K) 2 ~ 4 



I 
K 



(1) 



where L is the contour length, P is the persistence length, and K is the elastic modulus. 
We obtained L = 559.7, P = 25.2 and K = 230.2 by fitting in the region x > 525. The 
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fitting curve /ewlc, plotted in dashed line in Fig. [TJ matches well with the FEC in this 
region. The fitting values are consistent with the theoretical ones, L = (N m — l)bo = 562.1 
and K = 2/^5 = 220. It is known that P is a sum of two contributions; one is the 
intrinsic persistence length P originating from the chain bare stiffness, and the other is the 
electrostatic persistence length P e owing to the Coulomb repulsion between monomers on 
the chain. In this study, we calculated Pq by equating the bending energy Pbend to ksT/2 
for a chain segment of length P with the bending angle equal to 1 rad. Pbend is equal to 
(k2(b /P ) 2 + A; 4 (6 /P ) 4 ) x (Po/b ) because there are Po/b monomers on the chain segment 
and the bond angle 9 at each monomer is it — (bo/Po) in average. By solving this equation, 
we have Po = 12.2, which gives an estimation of P e equal to P — Pq = 13.0. It is worth 
noticing that Po obtained here is approximately equal to the major radius Rq of the toroidal 
condensate before stretching. This is a result of charge neutralization. The condensed 
tetravalent counterions neutralize the chain charge and therefore, the electrostatics does not 
play a major role in determination of the toroid size. At this moment, P e is roughly zero as 

n 

shown in the previous study of flexible chains [2J| . So the major radius relates directly the 
intrinsic persistence length of a PE. 

We calculated the excess work for the toroidal condensate AW = / 60 (/ — /ewlc)^ an d 
found AW = 1658. Therefore, the condensation energy is 3.47 per monomer (because the 
toroid part of the chain contains 478 monomers before stretching). This energy is mainly 
determined by the energy needed to break an electrostatic binding between a condensed 
tetravalent counterion and a monomer, which reads as |Ab&bX'(+4)(— 1)/1.5| = 14.98 or 
equivalently 3.74 per monomer. This finding suggests an electrostatic origin for the PE 
chain condensation. 

The condensed chain preserves the rod-toroid-rod structure in most of the time during 
stretching. In order to understand more deeply the elastic response, we calculated the 
average bond length on the toroid part, b t , and on the rod part, b r , of the chain. The results 
are presented in Fig. EJ^a). We found that bt stays roughly a constant up to x ~ 450, and the 
value is 1.08. This value is smaller than the equilibrium value bo = 1.1 of a bond, obviously 
due to the condensation of the tetravalent counterions which tightly bind middle part of the 
chain to form a toroidal structure. It is this tight binding which makes shorter the bond 
length b t . Since the toroid part consists of several loops, the effective elastic modulus on 
this part is very stiff. Therefore, the variation of b t with the extension x is hardly seen. On 
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FIG. 2: (a) bt and b r , (b) L c s as a function of x. W is plotted in the figure (a) to indicate the 
occurrence of the loop-by-loop unfolding. 



the other hand, b r on the rod part displays a sawtooth structure, coinciding with the FEC. 
This shows that the elastic response of the chain is mainly determined by the rod part. If / 
is strong enough to overcome the binding between the strands on the toroid, a loop is pulled 
out of it. 

We also calculated the effective contour length of the chain L e s, defined as the total 
length on the rod part plus the diameter of the toroid. The result is shown in Fig. [2(b). 
We see that L e g, on the main trend, depends linearly on x, but, locally, shows a series of 
small peaks (enclosed inside circles in the figure). Each peak corresponds to the moment 
when a loop is pulled out of the toroid and the tension of the chain is suddenly released; 
consequently, / is sharply decreased and so is b r . The linear dependence of L cS tells us that 
x/L e g is approximately constant during the pulling process, which suggests the existence 
of a reference force fo. fo can be estimated by replacing x/L in Eq. (CD) by x/L e g. Since 
x/L e ff ~ 0.95 in this study, fo is about 3.5. 

We did find that the sawtooth structure of FEC oscillates around this reference force. 
These oscillations can be attributed to the fluctuations of L cS (or b r ) due to the loop-by- 
loop unfolding of the chain. This phenomenon has been observed experimentally jj, 4, ||. 
The observed reference force can be used to estimate x/L e g in experiments by x/L c s = 
1 — \JksT ' I '4P fo jsj, which gives a typical value of 0.8 for DNA condensation, for e.g., by 
spermidine. This value is smaller than our simulation value 0.95. The discrepancy is due 
to the higher counterion valency and the smaller ion size used in this study, which leads 
to a stronger binding force between toroid loops than in the real experiments and hence, 
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a larger f . Consequently, the condensation energy obtained here is higher than a typical 
value obtained by experiments, 0.08 to OSksT per base pair 

We further investigated how the size of the toroid varies under stretching and the results 
are shown in Fig. [3j We observed that the major radius Rq decreases linearly until the 
release of the first loop. It abruptly increases at the moment when a loop is pulled out of 
the toroid, and then decreases with the extension. This behavior repeats many times and 
the final curve shows a sawtooth structure. For the minor radius r , we found that it is a 
downstairs function, decreasing simultaneously with W. When a loop is pulled out of the 
toroid, the number of the chain strands in the annular tube of the toroid decreases by one, 
and thus r exhibits a discontinuous jump. These findings show that at the first stage of 
stretching, the toroid decreases its size constantly to some value. The size then oscillates 
around this value (i?o — 5 in this study) while the toroid starts to loose its loops, one by 
one, during the stretching process. 

This structural transition has been approved by the snapshots, shown in Fig. H](a). The 
snapshots clearly show that the diameter of the toroid firstly decreases to a critical value 
(pictures (1) to (3)) and then keeps the value around this value for the rest of the loop 
releasing process (pictures (4) to (8)). The breaking of the toroid into several loops (or 
toroids) is not observed in the simulations. We have verified that the other independent 
runs show the consistent elastic response. The W curves for these 20 independent runs 
are plotted in Fig. 11(b). The consistency of these curves confirms that the whole process 
is repeatable and the nonequilibrium effect has been minimized by our slow pulling speed. 
Therefore, the results can be used to understand the elastic response in a near-equilibrium 
condition. We, furthermore, calculated the mean value of x at which W decreases by a step. 



9 



(1) W=6.5 




(7) W=2.5 



(8) W=1.5 



(9)W=0 



(b) 




100 



200 



300 



400 



500 



FIG. 4: (a) Snapshots of simulations. The toroid part of the chain is colored in light (green) color, 
(b) W as a function of x for the 20 independent runs. The inset shows the averaged position (x) 
to occur a loss of one loop under stretching. 



The results, W vs. (x), are plotted in the inset of Fig. S](b). We see that (x) lies on a straight 
line, which demonstrates that the circumference of a loop is a constant. The slope of the 
line corresponds to a circumference A (x) = 37.03 per loop, or equivalently a circle of radius 
5.9. This radius is consistent with the value of Ro at the peaks of the sawteeth in Fig. [31 

It is known that the persistence length of DNA is about 50nm. Consequently, when 
DNA is collapsed into a toroid, the periphery length of the toroid is approximately 300nm. 
Murayama et al. [5:] have shown that the stick-release pattern has a periodicity characterized 
by a distance of 300nm. This result supports the loop-by-loop unfolding of DNA, as do 
our simulations. Moreover, they found that the higher the salt concentration, the earlier 
the stick-release pattern happens in the stretching. Thus, the chain looses its first loop 
at a smaller extention position, according to our results. This can be understood by the 
phenomena of reentrant transition, in which a condensed DNA becomes less and less stable 
while the salt concentration is increased [§, [lsj. Because the loss of the first loop of toroid 
occurs at a larger R at higher salt concentration, we predict that the periodicity of the 
stick-release pattern increases with salt concentration. 
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FIG. 5: 9 and both in unit of angle degree (°), as a function of x. W is plotted to indicate the 
occurrence of the loop-by-loop unfolding. 



Finally, we studied the zenith angle 9, measured from the pulling- force direction to the 
normal direction of the toroid surface, to give insight of the unwrapping pathway of a toroidal 
condensate under tension. 9 is plotted in Fig. [5] as a function of x. We found that the chain 
was stretched from a starting zenith angle of about 65°. In the stretching process, 9 showed 
piecewise decreasing behavior with x, happened coherently with the change of the winding 
number. Therefore, the toroid normal was pulled, over again, toward the force direction 
during the stretching and swung back to its early direction at each moment when the toroid 
released a loop. It is the sudden drop of the chain tension which restores the normal 
vector back to its early direction. We also calculated the azimuthal angle <fi of the normal 
vector, plotted in Fig [5] too, which shows that the projection of the vector on the plane 
perpendicular to the pulling direction fluctuated around a reference direction. Kulic and 
Schiessel predicted a flipping-up-and-down pathway of transition for the normal vector of a 
spool of DNA helically coiled onto a histone protein {l^ . Such transition was not observed 
in this study. This finding suggests an important role played by the histone protein. For 
a DNA-histone spool under tension, the radius of the spool is restricted by the histone 
core. To pull a loop out of the spool inevitably induces a flipping-up-and-down movement 
of the spool. On the other hand, for a spool of DNA wrapped onto itself by multivalent 
counterions, a loop can be released gradually by decreasing the radius; therefore, the motion 
of the spool is more gentle and 9 oscillates in a small range of angle, as seen in our case. 
This delicate difference of the dynamics between these two systems will be clarified more 
clearly in the future. 
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In summary, we have demonstrated that the loop-by-loop unfolding of a toroidal PE 
produces the stick-release pattern in FEC. We have shown how the internal structure of a 
toroidal condensate and its normal vector change during stretching. The results give deep 
insight of the elastic response and the structural transition of condensed DNA molecules 
being stretched. 
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